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Abstract 

We determine the infinite volume coefficients of the perturbative expansions of the self-energies 
of static sources in the fundamental and adjoint representations in SU(3) gluodynamics to order 
a 20 in the strong coupling parameter a. We use numerical stochastic perturbation theory, where 
we employ a new second order integrator and twisted boundary conditions. The expansions are 
obtained in lattice regularization with the Wilson action and two different discretizations of the 
covariant time derivative within the Polyakov loop. Overall, we obtain four different perturbative 
series. For all of them the high order coefficients display the factorial growth predicted by the 
conjectured renormalon picture, based on the operator product expansion. This enables us to 
determine the normalization constants of the leading infrared renormalons of heavy quark and 
heavy gluino pole masses and to translate these into the modified minimal subtraction scheme 
(MS). We also estimate the four-loop /3-function coefficient of the lattice scheme. 

PACS numbers: 12.38.Gc,11.15.Bt,12.39.Hg,12.38.Cy,12.38.Bx 
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I. INTRODUCTION 



Perturbative expansions in the coupling parameter a of four dimensional non-Abelian 
gauge theories, Yl^=o c n aTl+1 y are expected to be asymptotic. The structure of the operator 
product expansion (OPE) determines one particular pattern of asymptotic divergence. This 



is usually named a renormalon |1| or, more specifically, an infrared renormalon. Its existence 
has not been proven. It could only be tested in QCD by assuming the dominance of /3o-terms, 
which amounts to an effective Abelianization of the theory, or in the two-dimensional O(N) 
model j^j], where it is suppressed by powers of 1/N. Moreover, the possible non-existence 
or irrelevance of renormalons in Quantum Chromodynamics has been suggested in several 
papers, see, e.g. Refs. (3), Q] and references therein. This has motivated dedicated high order 
perturbative expansions of the plaquette, see, e.g. Refs. p-|8|, in lattice regularization, with 
conflicting conclusions. Powers as high as a 20 were achieved in the most recent simulation [9] . 
However, the expected asymptotic behavior was not seen. A confirmation of this "non- 
observation" in the infinite volume limit would significantly affect phenomenological analyses 
of data from high energy physics experiments where renormalon physics plays a fundamental 
role. This is certainly so in heavy quark physics, where addressing the pole mass renormalon 
is compulsory for almost any precise computation, such as for determinations of the heavy 
quark masses in the MS scheme, the decay of heavy hadrons, or heavy quarkonium physics. 

Fortunately, in a recent letter, the existence of renormalons in quantum gluodynamics has 
been unambiguously established [lo|]. The quantities studied were the self-energies of static 
sources in the fundamental and adjoint representations. This analysis clearly identified the 
reasons for the previous non-detection of the renormalon-associated asymptotic behavior of 
the plaquette. In lattice regularization with the Wilson action, renormalon dominance only 
sets in at very high orders in n. In the case of the static self-energy, an operator of dimension 
d = 1, the renormalon behavior was confirmed for n > 9. Therefore, for the plaquette and 
the associated gluon condensate, an operator of dimension four, we expect n ~ 4 x 9 to 
be necessary to confirm the expected asymptotic behavior, an order that is quite beyond 
those reached so far in simulations. On top of this, it was shown that for presently reachable 
volumes the proper incorporation of the leading finite size effects (FSE) is required to obtain 
the correct infinite volume limit, something that had not been done previously either. Finally, 



in Ref . 



10[ preliminary results for the normalization of the renormalon were obtained, which 



turned out to be perfectly consistent with expectations from continuum computations in the 
MS scheme. In this article we provide greater detail on these simulations and our analysis 
methods, present finalized results, and further extend this previous study. 

This article is organized as follows. In Sec. [Ill we review numerical stochastic perturba- 
tion theory (NSPT), our improvements on previously existing techniques, and the specific 
aspects of the lattice computation relevant for our case. In Sec. II I II we define our primary 
observable: the self-energy of a static source, and detail the expected asymptotic behavior of 
its perturbative expansion due to the leading renormalon. In Sec. II VI we define the Polyakov 
loop, relate this to the static self-energy, and explain how our primary data sets are ob- 
tained. In Sec. |V]we present a theoretical study of the leading FSE and how these will affect 
the signatures of renormalon dominance. Subsequently, in Sec. [VT| we investigate, mostly 
numerically, subleading FSE that may pollute our data and estimate their systematics. In 
Sec. IVHI we determine the infinite volume coefficients of the perturbative expansion, study 
their renormalon structure, and extract universal results in the lattice and MS schemes, 
before we conclude. 



II. LATTICE IMPLEMENTATION 

Below we discuss the simulation method and its implementation. After a brief intro- 
duction into NSPT we detail a new second order integrator, introduce twisted boundary 
conditions and link smearing. 



A. Stochastic Quantization and NSPT 



Stochastic Quantization (SQ) llj enables the calculation of expectation values in quan- 
tum field theories and presents an alternative to, for instance, the path integral formalism. 
In recent years, SQ was employed in several studies within different fields of physics, ranging 
from the quark-gluon plasma 12|. even addressing the notorious sign problem of QCD at 
non- vanishing baryon densities \l<L \h\ , to quantum gravity jl5| . SQ turns out to be effi- 
cient also from the point of view of computer simulations due to the absence of any global 
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accept/reject step thus allowing, in principle, for a fast update of the system under consider- 
ation. The draw-back is the requirement to span a range of integration step sizes, to enable 
an extrapolation to continuous stochastic time. 

For simplicity, we assume a scalar field <p(x) depending on spacetime x and dynamics 
governed by an action S[<f)]. The core of SQ, the Langevin equation, then reads 

d<f ) (x,t) _ dS[<j>] 

at ~ d<f>(x,t) v[ , ) ' [ ) 

where t is the so-called stochastic time. The T](x,t) is a Gaussian noise variable with the 
properties 

{7](x,t)) v = 0, 

(r](x, t)r]{x\ t')) v = 25(x - x')5(t - t') . (2) 

The subscript "77" stands for an average over the noise. Given a generic observable A(4>), it 
can be showrj^ that the time average 

A$)= lim i fdtA^) (3) 

coincides with the expectation value on the quantum vacuum, i.e., 

1 



A(0) = - J[d<t>]A{<t>)e- s W, (4) 

where Z is the partition function. 

If the degrees of freedom of the system under consideration are not scalar but obey a 
group structure, as it is the case for lattice QCD, the above machinery has to be modi- 



fied accordingly (numerical stochastic perturbation theory, NSPT 17|, ll8[, for a review see 



Ref. |19|). In lattice simulations, spacetime is discretized by introducing a four- dimensional 
hyper-cube of iVj x Nt sites, where asymmetric volumes N$ 7^ Nt are legitimate. A pe- 
culiarity of NSPT is that no mass gap can be generated in perturbation theory. Hence the 
lattice spacing a is neither set nor determined a posteriori, so any NSPT-related reference to 
a is purely formal. For instance, the limit N$ — > 00 can either be interpreted as the infinite 
volume limit L = aNs — > 00 at fixed a or as the continuum limit a — > at fixed lattice extent 



1 For a proof in perturbation theory, see Ref. [16 1. 
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L in physical units. Lattice sites n are referenced by their spatial and temporal coordinates, 
Hi E {0, ... , N s — 1} and n 4 £ {0, . . . , N T — 1}, respectively 

The gauge degrees of freedom A^(x) in the continuum are elements of the Lie Algebra of 
SU(3) in representation^ R. On the lattice these are implemented as compact link variables 
U^(n) ~ e* A ?' (n+1 / 2 - )a l £ SU(3), connecting the sites n and n + ft, where fx denotes a unit 
vector in direction ll. 

The straightforward generalization of the Langevin equation Eq. ([1]) to fundamental link 
variables reads 

J^7 M (M) = -tY,T a [K,,S[U] + V ;(n,t)]U,(n,t) , (5) 

a 

where S[U] is the gauge action and T a , a = 1, ... ,8 are the traceless Hermitian generators 
of the SU(3) Lie algebra with the normalization ti{T a T b ) = \5 a b- We define the derivative 
within Eq. (JSJ) of a function f(U) with respect to a Lie group variable U following Ref. [20 ]: 

/ (e'Ea^U) = f(U) + J2co a V a f(U) + 0{u 2 ) , (6) 

a 

where u a are small real parameters. 

Perturbative lattice simulations up to n max loops become possible by a formal weak 
coupling expansion of the gauge fields. In the algebra and group this reads 

A = A (1) /H + A^/r 1 + • • • + A (2nmax) /T nmax , (7) 
U = 1 + + U^/3- 1 + ■■■ + u {2nm ^ f3~ n ™* . 

Above, (3 denotes the lattice coupling and relates to the strong coupling parameter as = 
g 2 /6 = (2vr/3)a. Note that while the A® belong to the Lie algebra of SU(3), the U {i) are 
no group elements. U however is, up to terms of 0(/H nm " an SU(3) group element. 
By Taylor expanding the exponent and logarithm of the two series, respectively, one can 
conveniently switch between algebra and group representations. Plugging the expansion 
Eq. (jTj) into a discretized version of the stochastic differential equation Eq. (j3J), one finds 
that the noise directly acts only on while the evolution of higher orders is governed 
by a hierarchical system of ordinary differential equations. In particular, the evolution of a 



2 Representation R has the dimension cIr. Here we consider two representations: the fundamental triplet 
(d,R = 3) and the adjoint octet (cZr = 8). 
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given order in stochastic time only depends on preceding orders 1, . . . ,i — 1 so that a 
truncation at finite 77, max is possible. 

The naive computational effort of NSPT scales like and the memory requirement 
like n max , compared to a factorial growth of the number of diagrams ~ n max \ in conventional 
perturbation theory. This makes high order expansions feasible. On an absolute scale, com- 
putation time of course becomes an issue for large lattice volumes or high n max , requiring 
optimizations of the NSPT algorithm. This study would have exceeded our present com- 
puter resources had we not used an improved numerical algorithm to evolve the Langevin 
equation Eq. (J5J). Its advantages were detailed in Ref. [2l|. Below we present the algorithm 
in detail. 



B. The second-order integration scheme 

The numerical integration of the Langevin equation Eq. (jSJ) requires the discretization 
of the stochastic time t, introducing a time step e (t = t m = me with integer m) and a 
prescription for the ^-derivative in Eq. (jSJ). Revisiting the scalar example, schematically the 
updating step for the zth degree of freedom 0, readqf 

0S m+1) = 0l m) - ft ] , (8) 

where the bracketed superscript labels the evolution in Langevin time t = me and /$ is a 
force term. In the simplest (Euler) integration scheme the force is given by 

// m > = e^S^ + V~evt } , (9) 

with the functional derivative V defined in Eq. (jBJ) for gauge theories and vff 1 — y/er)(n, t = 
me). 

Information on how the discretization changes the equilibrium distribution relative to the 
continuous-time expression of Eq. (SJ) can be drawn from the Fokker-Planck equation. To 
work this out, we label the probability distribution after m + 1 updates as V^ m+1 \(p): by 
defining W{<f>' 0) as the probability of jumping from configuration <fi to configuration 0', 

3 We have replaced the dependence on discretized spacetime coordinates n by an index i for simplicity. 
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we obtain the equality 

V (m+1 \<f)') = l[d<P} W{<j>' <- 0)P (m) (0) 



= (wm n w - 0, + /,)^ (m) (0) • (io) 

i 

The above product extends over all degrees of freedom and we have rewritten the probability 
of moving from <p to <$' in terms of ^-functions, involving the noise (that is implicit in /j). 
After some algebraj^ one obtains 

oo 1 

v^ +l \4>) = v^\ct>) + -Vn • • • v b . [<& . . . 4>„P(0)] . (11) 

3=1 J ' 

We recall that p( m+1 )(0) = V^ m \cj)) at equilibrium, insert force terms into Eq. (fTTj) and 
expand with respect to e. This leads to the identity 

= V,{V^[0] + V 4 }P(0) , (12) 

whose solution reads V(<j>) oc e~ s ^ with 

5[0] = S[0] + e5 1 [0] + e 2 5 2 [0] + .... (13) 

Within the above equation, S[<f>] is the original action of Eq. ([TJ). Thus, the correct equilib- 
rium distribution and, consequently, Eq. (J5J) is recovered in the limit e — >• 0. In the Euler 
scheme, for example, S[(f)) is given by 

S[<j>] = S[<f>] + € -J2 (WiViS^] - ViS[<l>}ViS[<f>)) + °^ 2 ) • (14) 

i 

We detailed the formalism for a scalar field <fi. In the case of non-Abelian SU(iV c ) gauge 
theory, the discretized Langevin update readjfl 

{/ KD =e - l OT [/ W ) (15) 

where the force term in the Euler scheme is given by the analogue of Eq. ([9]): 

ft = e^S\U] + y/irff. (16) 



4 Essentially, one represents each ^-function as a Fourier integral, Taylor-expands in the force term, expresses 

each power of the expansion by suitable derivatives with respect to the </>s and integrates by parts. 

5 The index i now contains both spacetime position n and direction /x. 
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With the group derivative defined as in Eq. fl6]) the above procedure can be repeated for non- 
Abelian degrees of freedom, again leading to a Fokker-Planck equation. The only difference 
lies in the fact that group derivatives do not commute. More precisely, in the continuum 



[v^,v^] = -r bc v c M,- 



(17) 



where f abc are the structure constants of the Lie algebra. Obviously, this has a non-trivial 
impact on the equilibrium distribution S[U] at e > 0. For instance, plugging Eq. (1TB]) into 
the Fokker-Planck equation, we obtain 



S[U] 



1 + 



^f) S[U] + Z E ( 2V ^S[U) - V*S[V\S7?SP\) + 0(e 2 ) , (18) 



where Ca = N c is the quadratic Casimir invariant of the adjoint representation of SU(iV c ). 

From Eqs. (fT4"l) and ffl8|) it is evident that numerical simulations with different values of 
e are necessary to extrapolate to continuous stochastic time e — > and to recover Eq. (jSj) 
and the continuum distribution. Simulations at small e obviously are more costly and it 
is tempting to keep e as large as possible. However, for large time steps corrections to the 
leading linear dependence will become sizable and extrapolations to e = less controlled. 

A reduction in computer time while maintaining a safe e — > extrapolation becomes 



possible by employing higher-order integration schemes. 



schemes exist up to the third order for Abelian theories 22 



b our knowledge, Runge-Kutta 



23[ (the general solution to the 



Fokker-P 
theories 



anck equation is known), and up to the second order in e for non- Abelian SU(iV c 



24 



25|. In the latter case, only one variant of the general solution is published, 



namely the two-step algorithm 



U 



U> = e -iJ2a T i^S[U]+V-zrit) f/M 5 
(m+1) = ~iJ2 a T a (^S[U]+^S[U']+^^S[V'}+^vt) rr(m) 



(19) 

(20) 



where S[U] and S[U'] stand for the action computed using the fie 



We refer to this second-order integrator as the "BF scheme" 



ds and U', respectively. 



251 ] . Note that the evolution 



cannot be factorized into sweeps involving single link updates: both U[ and have to 

be computed for all links i, prior to the replacement of the original field u\ m \ In particular, 
in the second step both S[U] and S[U'] are needed. This requires three copies to be kept in 
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memory concurrently of 2n max + 1 orders of complex three by three matrices for each lattice 
link. 

Below we derive the general solution and provide an optimized alternative to Eqs. (fT9l) 
and (1201) which not only saves matrix additions but also reduces the memory requirements. 
The general ansatz for the second-order algorithm reads 

JJ> = e iE a T a (kieVfSlU]+k2V~eV?) Jj( m ) ; (21) 
jj(m+l) = e -iX; a ^(fc3eV i a S[C/]+fc4eV i a S[U']+A :6 C yl ^V < a S[[/']+fc6V^<) Jj( m ) (22) 

Plugging the force term of Eq. ( 122|) into the Fokker-Planck equation and Taylor-expanding 
the derivative of S[U'], after some algebra some constraints are obtained: at O(e ) the 
non-Abelian analogue of Eq. ( ITgj) yields 

k 3 = l-k 4 , k 2 e = l, (23) 

in order to recover the correct e — > distribution, while the elimination of terms proportional 
to e (using k 3 = 1 — &; 4 ) results in 

l-4fc 4 ±2 v /2£; 4 (2£;4-l) fOA , 
kl ~ 2k d ' [M) 



h ~ 2k~k & ' (25) 

^ = -i ±M 3^3iy [ (26) 

where fc 4 and fee = ±1 can be chosen freely. The BF scheme is recovered setting fc 4 = | = k 3 , 
he = 1. The choice k 4 = k$ = I, however, further simplifies the algorithm: 

^(m+1) = ~i !Z a T a (e^S[U']+ ^- 3 ^ iV^Sm+^nt) ^(m) (2g) 

The gain of this variant is twofold: besides saving a matrix addition when computing the 
force term of Eq. f l2"g|) instead of Eq. fl2"U|) . there is no need to store (or to recompute) WiS[U]. 
After the intermediate step the original fields (and/or Si [77]) can be overwritten, reducing 
the memory requirement by one third. 

We tested the integrator defined through Eqs. fl27|) and f )28|) within NSPT: due to the 



need to rescale the time step e h->- e//3 (for details see, e.g., Ref. [19]), after inserting the 
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TABLE I. Comparison on N 4 lattices between analytical and NSPT results for the one-loop coef- 
ficient pq of the plaquette: The e-values used in the extrapolations range from 0.04 to 0.07. The 
analytical result is given by 2(1 — iV~ 4 ). 



N 


analytical 


Euler 


2nd-order BF 


new 2nd-order 


4 


1.9921875 


1.9931 (6) 


1.9922 (9) 


1.9924 (7) 


6 


1.9984568 


1.9985 (3) 


1.9987 (3) 


1 9986 C\) 


8 


1.9995117 


1.9997 (2) 


1.9992 (3) 


1.9996 (3) 


10 


1.9998000 


1.9996 (1) 


2.0001 (2) 


2.0001 (2) 


12 


1.9999035 


1.9998 (1) 


1.9999 (1) 


1.9998 (1) 


TABLE II. Comparison on N 4 lattices between diagrammatic lattice perturbation 


theory (DLPT) 


and NSPT results for the 


two-loop coefficient p\ 


of the plaquette. 




N 


DLPT 


Euler 


2nd-order BF 


new 2nd-order 


4 


1.20370366 


1.2020(15) 


1.2005(17) 


1.2012(17) 


6 


1.21730787 


1.2173 (7) 


1.2166 (8) 


1.2180 (8) 


8 


1.21965482 


1.2203 (4) 


1.2178 (7) 


1.2199 (6) 


10 


1.22031414 


1.2203 (3) 


1.2204 (4) 


1.2212 (6) 


12 


1.22055751 


1.2208 (2) 


1.2200 (3) 


1.2204 (2) 



perturbative expansion into the discretized Langevin equation, it turns out that the contri- 
bution proportional to Ca in the force term of Eq. ([28]) only affects the two-loop level and 
beyond. 

In Tables HI and HT1 we compare one- and two-loop plaquette coefficients po and p\, defined 
through 

([/□) = l-Po/r'-Pa/r 2 ---- . (29) 

These were computed using the new second-order algorithm, diagrammatic lattice perturba- 
tion theory, the Euler integrator and the BF scheme for different symmetric volumes iV 4 with 
periodic boundary conditions (PBC). For the Euler integrator, the fit function employed in 
the extrapolation is constant plus linear while for the second-order schemes the ansatz is 
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constant plus quadratic (in the latter case, we checked that the coefficients of terms linear 
in e indeed vanish within errors when using a linear plus quadratic fit function). In all the 
cases we find agreement between the methods within two standard deviations. 



O 
Q_ 




0.002 0.004 0.006 0.008 0.010 

2 



FIG. 1. The one-loop plaquette coefficient p$ vs. e 2 for N 4 = 10 4 : Euler integrator, BF scheme and 
the new second order integrator. 



Figures [T] and [2] illustrate the finite-e plaquette results for the three integrators: while 
all sets extrapolate to the same limit within error bars, the ones corresponding to the 
new second-order scheme are clearly much flatter than the others. In particular, the e 2 - 
dependence of the new second-order integrator is greatly reduced compared to the BF scheme 
of Eqs. ffl9|) and (120]) . Note that we allowed for cubic terms in the curves drawn for the sec- 
ond order integrators. We will see in Sec. [IV] below that for the observables of interest in 
this work extrapolations in e 2 are so flat that in most cases a non-trivial slope cannot be 
resolved within statistical errors. 
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0.002 0.004 0.006 0.008 0.010 
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FIG. 2. The two-loop plaquette coefficient p\ vs. e 2 for N 4 = 10 4 : Euler integrator, BF scheme 
and the new second order integrator. 

C. Twisted Boundary Conditions 



Instead of PBC, one can also impose twisted boundary conditions (TBC) |26M29l] 



U^n + N s 9) = tt u U^n)ttl 
U,{n - N s 0) = QtU,(n)Q u 



(30) 
(31) 



where the links that pierce a lattice boundary of a twisted (spatial) direction v are multiplied 
by so-called twist matrices Q u that must satisfy 



n 3 , 



-if-H. 



(32) 
(33) 



Here z E {1, e i27r / 3 l, e i47r/3 l} is an element of the center of SU(3). The condition Eq. (|52|) 
guarantees that the value of the transported link (n + Ngfl + NgO) is independent of the 
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order with which two twisted boundaries /i, v are transversed. Gauge transformations A(n), 
which rotate the link variables according to 

U^n) ^ A(n)t^(n)A t (n + p) (34) 

must obey the same TBC Eq. ( J30l) . 

The measure as well as Wilson loops without net winding numbers across boundaries 
(such as the elementary plaquette within the action) are invariant under the transformation 



Ufj,(n) — > zUfj,{n) , Wn G {n : n ■ ft = const.} . 



(35) 



TBC rely on this center symmetry of the SU(3) gauge action and measure, and can be 
implemented for the link update either by multiplying the plaquettes in corners of twisted 
hyper-planes with suitable center elements, or by imposing Eq. ( l30l) with an explicit choice 
of Q u . We implemented the latter using 



/ n 1 n \ 



fix 



1 

1 

1 



' C* o 



/ o o c* ^ 



^3 — 1^2^i 



1 



(36) 



1 

V o ocy 

where £ = e 2l7r//3 , (* = l/£. This choice is arbitrary up to global unitary transformations: as 
long as Eq. f[3"21) is satisfied, the resulting physical amplitudes will not depend on the explicit 
choice of Q u . As the subscripts indicate, we impose the twist for all spatial directions. 
Twists in two directions have a non-trivial effect too, while twisting only one direction can 
be absorbed into a re-definition of the link variables. The effect of twist is twofold: TBC 
eliminate zero modes which otherwise require an explicit subtraction [3]. Furthermore, at 
least at low orders in perturbation theory, TBC reduce finite size effects as the possible gluon 
momenta are restricted to integer multiples 28j| of 



Pu 



, v = twisted direction . 

2tt 



(37) 



x . ^ ■ = periodic direction . 

This means that gluon momenta in twisted spatial directions reach values as low as 27r/ (3Ng), 
compared to 2tc/Ns in periodic directions. So, roughly speaking, the modes in a twisted 
direction behave as if the corresponding lattice extent was 3N$ instead of N$. We refer to 
the cases of twists applied to two and three directions as TBCxy and TBCxyz, respectively. 
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TABLE III. Plaquette coefficients. 4 4 PBC: DLPT (first two orders) and NSPT (remaining or- 



ders iMj)- 4 4 TBCxy and TBCxyz: DLPT (first 



order) and NSPT (remaining orders). 32 4 PBC: 



DLPT (first two orders) and NSPT (O^' 6 ) fig]). Infinite volume: DLPT from Ref. jMj, using 



the lattice integrals of Ref. \32l 1. For all NSPT data, the e — > extrapolation was carried out. 



order 


4 4 PBC 


4 4 TBCxy 


4 4 TBCxyz 


32 4 PBC 


oo 4 


r 1 


1.9921875 


2 


2 


1.9999981 


2 


/r 2 


1.2037037 


1.2184(5) 


1.2200(3) 


1.2207904 


1.2279575 




2.887(3) 


2.955(2) 


2.957(2) 


2.957(3) 


2.9605(1) 


r 4 


-9.05(1) 


-9.41(1) 


-9.40(1) 






r 5 


-32.49(6) 


-34.51(9) 


-34.34(5) 







The effect of TBC is noticeable in particular on small lattice volumes, as Table II III il- 
lustrates for the average plaquette. The two- and three-loop TBCxy and TBCxyz data 
obtained on 4 4 volumes are close to the infinite volume (as well as to 32 4 PBC) results at 
two and three loops. This clearly is not the case for 4 4 PBC data. Note that both analytical 
one-loop TBC coefficients happen to be volume-independent on symmetric lattices, due to 
cancellations between different plaquette orientations. 

The situation is different for the Polyakov loop L defined in Eq. (J63l) below. First of all, 
for this observable it matters whether it is obtained in an untwisted or a twisted direction. 
We calculate L in untwisted directions, for which no modification is necessary with respect 
to PBC, and extract the static energy 5m via Eqs. ( 164^) - (!65|) . As it was shown in Ref. 33 j 
for this observable, TBC significantly reduce FSE, resulting in a much flatter extrapola- 
tion towards infinite volume. If this flatness at low orders was taken as the only criterion, 
TBCxy would be the boundary condition of choice. However, it turns out that TBCxy has 
a draw-back compared to TBCxyz: in non-perturbative simulations only the latter prevents 
tunnelin g be tween different Z(3) phases while TBCxy merely leads to a reduction compared 



to PBC 34]. As a consequence, small volume TBCxy simulations were found to fluctuate 
more and to return noisier signals than their TBCxyz counterparts. Regarding the statistical 
fluctuations we made a similar observation, even though tunneling between Z(3) sectors is 
not an issue in our NSPT simulations since = 1 . Fig. |3] shows stochastic time histories 
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FIG. 3. Stochastic time series of one-loop (Lq, left) and 12-loop (L\\, right) Polyakov loop coeffi- 
cients for TBCxyz (blue) and TBCxy (red) on 16 lattices for e = 0.05. 

obtained on 16 4 volumes at fixed e = 0.05 for TBCxy and TBCxyz of the one-loop and 12- 
loop coefficients of the Polyakov loop. While the trajectories of the one-loop coefficient L 
show a similar behavior for TBCxy as for TBCxyz, we observe a peak in the twelve loop L\\ 
TBCxy measurement history, which is symptomatic for TBCxy simulations. The enhanced 
numerical stability and smaller fluctuations, in particular at large orders of expansions, mo- 
tivated us to choose TBCxyz for this work. A better understanding of the origin of these 
differences between TBCxy and TBCxyz would be desirable. 



D. Link smearing 



The lattice discretization of observables and action is not unique. For instance one can 
construct Wilson loops and Polyakov loops, replacing the link variables U^(n) by fat or 
"smeared" links. In the context of a lattice determination of static potentials and of static- 
light meson masses this was for instance done in Ref . 35j , to reduce the self-energy, enabling 
an improved signal to noise ratio at large Euclidean times. As long as the smearing is an ultra- 
local procedure, defined on the scale of a few lattice spacings, making this replacement in a 
Polyakov loop corresponds to a different choice of discretization of the static action. Smearing 
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is sometimes also used within the definition of fermionic actions, see, e.g., Refs. 

Several smearing methods are available, one of which is known as analytic or "stout" 
Stout links are automatically elements of the SU(3) group, without a numer- 



smearmg 
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J2j,(n) - R,(n) - ~tr (4 ( n ) - 1 



(39) 



ically delicate projection into the group. Therefore, implementing stout smearing within a 
perturbative expansion is straightforward. Stout smeared links are obtained by the replace- 
ment 

UM ^ U s ; out (n) = exp(zQ»)^(n), (38) 

where Q^{n) is Hermitian and traceless and hence in the algebra by design: 

i 

i^(n) = C»f/» , (40) 
C» = Yl P^UMU^n + £>)Z7t(n + p) + Ul(n- 0)U,(n - u)U u (n - u + /})) . (41) 

Note that within the sum of staples C^(n), surrounding the link U^n), the sum convention 
is not implied and p^ u are weights that can be set at will. In our case, we choose pi V = 
and pa = p = 1 /6 otherwise. The value of the weight was chosen to minimize the one-loop 
static self-energy after one smearing iteration. We remark that this is not necessarily the 
best possible choice, e.g., in a non-perturbative setting. We apply only one smearing step to 
keep the static action local. 



III. SELF-ENERGY OF A STATIC SOURCE 



In this section we introduce our conventions, relate self-energies of static sources to heavy 
quark and heavy gluino pole masses, and discuss the expected renormalon structure. 

The triplet and octet self-energies are defined as the lowest energy eigenvalues of the ef- 
fective Hamilton operator in temporal gauge of the sector of Hilbert space of gauge triplet 
and octet states with respect to gauge transformations, applied to a fixed position. It is 
impossible to obtain the continuum limit for these self-energies (irrespectively of the repre- 
sentation), as these will diverge linearly with the ultraviolet cut-off |] Therefore, the value of 
the static self-energy depends on the chosen regulator. Yet, any hard-cut-off regularization 
scheme is suitable for the following discussion. In this article we use lattice regularization 
and write the self-energies in the fundamental and adjoint representation in the following 

6 In dimensional regularization this object is exactly zero, since the ultraviolet and infrared divergences 
(infrared and ultraviolet renormalons) are regulated by the same factorization scale and their sum vanishes. 
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way: 

1 oo 1 oo 

Sm = -J2 c^ p) a n+l {l/a) (fundamental), 5m~ g = - ^ c^ p) a n+1 {I / a) (adjoint) , (42) 

n=0 n=0 

where a -1 , the inverse lattice spacing, provides the ultraviolet cut-off. The coefficients c n 
depend on the details of the regularization, i.e., on the action used and the specific definition 
of the static propagator. We only consider the Wilson action 40[ but we explore two different 
definitions of the static propagator, with smeared (p = 1/6) and with the original (p = 0) 
temporal links. We label this dependence with a generic p G {0, 1/6}, see Sec. Ill Dl above. 

One may suspect that the dependence on the regulator might turn this object uninter- 
esting from the theoretical point of view. This is actually not the case since, for large n, c n 
becomes regulator independent, universal and equal to r n /u, the n + 1 order coefficient of 
the perturbative expansion of the pole mass 



^os 



(43) 



n=0 



up to (9[exp(— l/n)]-terms (due to subleading renormalons) . On an intuitive level this is 
clear, as the static energy and pole mass share exactly the same infrared behavior (up to 
0(l/m) corrections), which should cancel in the difference. 

The asymptotic behavior of r n can be determined assuming that the perturbative series is 
asymptotic and the validity of the OPE. The running of a{y) is governed by the /3- function 



(3(a) 



da 



where in our normalization 



-2a 



a \ 



( a \ 3 



\AttJ 



+ ... 







11 



ft = 102 , /?. 



MS 



2857 



0: 



latt 
2 



-6299.8999(6) . 



(44) 



(45) 



From 02 onwards the coefficients depend on the scheme. While is known [41], in the 



lattice scheme only (3^ u has been computed 32 
constants 

, Pi 



42 



43]. For convenience we define the 



si 



S-2 



2/r 

PI - A)& 
46/3 4 ' 

g + Affifafo - 2foPlP 2 + P 2 (-2[3f + g) - 2/3 4 /3 3 
326(6 - l)/3$ 



(46) 
(47) 
(48) 
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where only b is scheme- independent. In a given scheme, the A parameter is defined as 

A= lim z, e -[ 61n (&)+J" aM ^)] (49) 

v— >oo 

2 



A)0(^) X 
47T 



27T V 27T 



and we use Aqcd ~ A^jg synonymously with the size of a typical non-perturbative binding 
energy. 

A simple scheme- and scale-independent observable is the B meson mass. In the heavy 
quark limit this can be decomposed into the 6-quark pole mass mos and the remaining energy 
from the light quark and gluon dynamics A B = CbAq CD : 

m B = m s + A B + ( j . (50) 

This relates to the fundamental representation. For the adjoint representation we can think 
of a heavy gluino attached to gluons: 



m d = m~ gjOS + A H + 0[ , (51) 

V m g,OS/ 

where Ah denotes the dynamical contribution of the gluons (and sea quarks) to the gluino- 
nium mass. For both representations the uncertainty of the perturbative series of the pole 
mass will be of O (Aqcd), the next term in the OPE, since only the sum of the pole mass 
and the binding energy has a physical meaning. This ambiguity results in the successive 
contributions r n a n+1 to decrease for small orders n down to a minimum at no ~ l/(|a<i|a), 
where = f3 /(2ird) with d = 1. After this order the series starts to diverge, so one neglects 
the higher order contributions and estimates the error by the size of this minimum term, 
r n a?10+1 ~ exp[— l/(a d a)} ~ Aq CD /m. If the perturbative expansion has an ambiguity of 
order Aq CD then d = n. To quantify this behavior it is convenient to consider the Borel 
transform of the above perturbative series 

oo 

dte- t/as B[m os }{t), B[m os ](t) = r n~ ■ (52) 

n=0 U - 







The behavior of the expansion Eq. fj43|) at large orders is dictated by the closest singularity 
to the origin of its Borel transform, which, for the pole mass, is located at t — 2n/Po, i.e., 
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at u = 1/2, defining u = ^p. More precisely, the behavior of the Borel transform near the 
closest singularity at the origin reads 



1 



B[m $](t(u)) = N m v (l _ +b (1 + si(l - 2u) + s 2 (l - 2u) 2 + ••■) + (analytic term), 



(53) 



where by analytic term we mean contributions that are expected to be analytic up to the 
next renormalon (u — 1). This dictates the behavior of the perturbative expansion at large 
orders to be 



n^oo AT ,'PoY I> + l + 6) ( 6 6(6 -I! 

''" = "' V 2n I ~" ' " — 7 — r_rT ' Sl 



r(i + 6) 



(n + 6) (n + 6)(n + 6-l) 



s 2 



(54) 



This expression can be obtained from the procedure employed in Ref. 44j. The si-term was 



44J , and the S2-term in Refs. 



45 



46]. 



computed in Ref. 

As we mentioned, the large-n behavior of Cn' p ^ is the same as that of r n up to 0(e~ 1 ^ n )- 
terms (due to subleading renormalons) . Therefore, using the same scheme for the expansion 
parameter a, we obtain 



z (3,p) n^oo N 



f3 y i> + i + 6) 



-Si 



6(6-1) 



-s 2 + - 



(55) 



2n J r(l + 6) V' (n + 6) (n + 6)(n + 6-l) 
Note that all the dependence on the regularization details (and in particular on p) vanishes. 
The normalization constant N m also determines the strength of the renormalon of the singlet 
static potential, through the relation 



2N m + N v 







(56) 



since these contributions cancel from the energy E(r) = 2m + V s (r) 47H49]. 
For adjoint sources we have 



c vV) Nm . ( S° 
9 \ 2n 



T(n + 1 + 6) 



1 + 



si + 



6(6 - 1) 



-s 2 + 



■ (57) 



T(l + 6) V (n + 6) (n + 6)(n + 6-l) 
Again, the dependence on the regularization details (e.g. on p) vanishes, however, the octet 
normalization is different: N m . ^ N m . Eq. ( 157)) corresponds to the renormalon of gluelump 
masses (actually N m , = —N\, where A^a is the strength of the gluelump renormalon asso- 
ciated to Ah) and can be related to the pole mass and adjoint static potential renormalons 
through the relation 

2N m + N Vo + N A = 0, (58) 
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since the energy E(r) = 2m + V Q (r) + Ah is renormalon free 50]. 

To eliminate the unknown normalization constants, we may consider ratios. In a strict 
1/n expansion we havdl 



(3,p) n 



•'n— 1 



(59) 



2tt 



1 + --^ 



r? 



+ -L[^ + fe ( fe -i)( Sl - 2s2 )]+o^ 



This expression holds in any representation. It is also independent of the renormalization 
scheme used for a. Keep in mind though that the explicit expression does depend on the 



scheme, starting from 02- Here we will mainly use ai at t, where S2 is unknown and 02 = 02 i 
but we will also consider the behavior of the perturbative series in the MS scheme. 

Assuming that the coefficients are dominated by the renormalon behavior, we can deter- 
mine the order Uq + 1 that corresponds to the minimal term within the pole mass perturbative 
series. Minimizing r n a n+l results in 



(n + b) 



A)Q! 

2tt 



exp 



2 (n + b) 



+ o 



(no + by 



This yields the minimal term 



7T 



r(i + b) 




N m A [1 + 0(a)] 



(60) 



(61) 



While it is evident to most readers, we wish to emphasize that the perturbative series 
that defines the pole mass cannot be resummed (not even in a Borel way). Therefore, it does 
not exist in a mathematical sense and no rigorous numerical value or error can be assigned 
to this object. The most one could do is to define a pole mass to a given (finite) order N + l, 



m, 



(TV) 
OS 



= Sr&o r ™ a ™ +1 ' which w iH then depend on N. By taking iV ~ tiq we minimize this 
dependence^! One can then estimate the uncertainty of the sum to be (see, for instance, the 



discussion in Ref. 



5l|) 



no \r no \a 



no+l 



23/2-6^3/2 ny | A 



(62) 



r(i + 6) A) 

Note that this object is scheme and scale independent (to the 1/n-precision that we employed 
in the derivation) because, even though the normalizations N m and N m . depend on the 
scheme, the products iV m A and N m ~A are scheme- independent. 



io|. 



7 This equation corrects a mistake in Ref. 

8 In practice one would round N = int(no), giving a slightly different value. 
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IV. THE POLYAKOV LOOP 



We obtain the coefficients c n from the temporal Polyakov loop on hyper-cubic lattices. 
We investigate volumes of Nt lattice points in the time direction and spatial extents of N$ 
points. We choose PBC in time and TBC in all spatial directions (TBCxyz, see Sec. Ill CI) , 
eliminating zero modes and improving the numerical stability. For test purposes, we have 
performed additional simulations with PBC in all spatial directions to 0(a 32 ) for a 4 3 x 8 
volume and to lower orders for the specific volumes listed in the second row of Table IIVI 



TABLE IV. The Ns(Nt) values of the PBC and TBC runs. The different geometries are grouped 
in terms of the orders of the respective expansions C(a™ max+1 ). 





0(a 3 ) 


0(a 4 ) 


0(a 12 ) 


O(o 20 ) 


PBC 




4(4) 




8 (8,10,12,14) 




5 (5,6,7,8,10) 


4 (5,6,7,8,10,12,16,20,24) 


6 (6,8,10,12,16) 


7 (7,8) 


TBC 






8 (12,16) 
10 (8,12,16,20) 


8 (8,10), 9 (12) 
10 (10), 11 (16) 






12 (16,20) 


16 (12,16,20) 


12 (12), 14(14) 



The Polyakov loop is defined as 

l( r \n s ,n t ) 



ni 2? d R tT 



'JVt-1 



n,4=0 



R, 



n 



(63) 



where 



n 



iA«[(n+l/2)a] 



e *» 



G SU(3) denotes a gauge link in representation R, connecting 



the sites n and n + //, G {0, . . . , N s — 1}, n 4 e {0, . . . , Nt — 1}. We implement triplet 
and octet representations R of dimensions d,R = 3 and 8. The link U^(n) appears within 
the covariant derivative of the static action x/jD^ip, acting in the following way on a scalar 
lattice field f{n): D A f(n) = [U 4 {n)f(n + 4) - U\{n - l)f(n - 4)]/(2a). This discretization 
is not unique and we may substitute U4 by another gauge covariant connection. We use 
singly stout-smeared 39( covariant transporters (smearing parameter p = 1/6, see Sec. Ill DT) 
instead of U^n) as a second, alternative choice, to verify the universality of our findings. 
We perturbatively expand the logarithm of the Polyakov loop 

\n(L^p\N s ,N T )) __ 

n=0 



p( R >p)(N s ,N 7 



aN 7 



J2c ( n R ' p) (N s ,N T )a 



n+l 



(64) 
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in order to obtain the static triplet and octet self-energies in the infinite volume limit: 

5m = lim P (3 ' p) (N s , N T ) , 8m~ g = lim P^ 8 ' p \Ns, Nt) , (65) 

Ns,N T ^oo Ns,N T ^oo 



c^= lim ci R >»\N s ,N T ). (66) 

Ns,Nt— >oo 



where (see Eq. (Hz 



The primary objects that we compute are the coefficients c n (Ns, Nt)- The sets of 
c n (Ns, Nt) obtained on different geometries are statistically independent of one another. 
However, for a given volume, different orders n will be correlated. In computations of ratios 
c n{Ns, N T )/c n _i(N s , N T ), as well as in fits, we take these correlations into account. NSPT 
enables us to calculate the coefficients directly, i.e., that neither the lattice spacing nor the 
strong coupling parameter a enter the simulation. We have realized a large variety of TBC 
geometries, listed in Table HVl in addition to the PBC test runs. Each coefficient c n depends 
on Ns and Nt but also on the time step e of the Langevin evolution (see below). In this 



paper we employ the variant of the Langevin algorithm introduced in Ref. [21( and explained 
in Sec JII Bl which only quadratically depends on e. The time series were analyzed following 
Ref. [521 , allowing us to process either single runs or to evaluate sets of "farmed out" Monte 
Carlo branches. Special care was taken to ensure that every individual history for each order 
was sufficiently long relative to the respective autocorrelation time to guarantee a safe error 
analysis. Branches that failed this test were removed from the data analysis. 

Very high statistics runs were performed up to 0(a 3 ) and 0(a 4 ), to check if the coefficients 
of logs extracted from the data are in agreement with our theoretical expectations, to detect 
signs of ultrasoft \ih(Nt / Ns)-terms (see Sees. IVl and IVT1 below) . and to compare with results 
from diagrammatic perturbation theory. 

The bulk of data are obtained up to 0(a 12 ) and 0(a 20 ). We have kept the Langevin 
time between two successive measurements fixed, adjusting n up d ~ 56/e where n up d is the 
number of updates performed in between two measurements. For the a 12 and a 20 runs 
between 35000 and 80000 measurements were taken, corresponding (for e = 0.050) to 4 x 10 7 
- 9 x 10 7 updates. The integrated autocorrelation times increase with the order of the 
expansion and with the lattice volume. For instance, the integrated autocorrelation time of 
Co varied from 2.4 (6 4 ) to 15 (16 4 ), while that of cn from 9 to 30, in units of n up d- For our 
highest order coefficient Cig we found the values Tint ~ 18 and 29 for ll 3 x 16 and 12 4 lattices 
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respectively. Since the e = 0.050 measurements are separated by 1120 Langevin updates, our 
largest r int = 30 value corresponds to more than 33000 such updates, still leaving us with a 
few hundred effectively statistically independent measurements. 

Quadratic extrapolations in the Langevin time step to e = on a 16 4 lattice (based on 6 e- 
values) and on a 6 4 lattice (4 e- values) were performed up to 0(a 12 ). The 16 4 extrapolation, 
normalized to the value obtained at e = 0.05 (10 3 e 2 = 2.5) is shown for the unsmeared 
triplet coefficients ch' ^ (16, 16) in Fig. HI up to n = 11. A similar picture arises for the 
6 4 volume: within two standard deviations the extrapolated values are all found to agree 
with the results obtained at e = 0.05. We remark that for both volumes the fit functions are 
quite flat, with very small slopes in e 2 . The same was observed for the smeared and the octet 
coefficients. Based on this experience, the remaining volumes are only simulated for e = 0.05. 
However, the time step scaling was only tested within certain errors that are similar in size 
as (and in some cases larger than) the statistical errors obtained for the various geometries. 
We obtain a relative systematic error for each order by adding the statistical error of the 
extrapolation and the difference of the extrapolated value from unity in quadrature. For the 
other geometries this (multiplied by the e = 0.05 coefficient) is then added in quadrature to 
the respective statistical error. For orders larger than n + 1 = 12, we linearly extrapolate 
in n the systematics found for orders n + 1 < 12, to obtain an estimate. This procedure is 
performed not only for the coefficients themselves but also for ratios of coefficients. While 
the systematic time step error is assumed to be uncorrelated, we keep track of the correlation 
between the statistical part of the errors for different orders n, obtained on the same volume. 



V. FINITE SIZE EFFECTS FOR N T -> oo 

The finite size effects of P(N$, Nj>) [Eq. 064jl ] are well suited to a theoretical analysis in 
the limit Nj- — > oo (actually > Ng for most of the geometries that we simulate). In this 
limit the self-energy of a static source in a finite spatial volume is obtained^ 

5m(N s ) = Jim P(N S ,N T ) and c n (N s ) = Jim c n (N s ,N T ) . (67) 



9 The discussion in this section applies to any representation and to smeared and unsmeared Polyakov loops. 
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FIG. 4. Time step extrapolations of ratios Cn' ^ (16, 16; e)/cn 3 '°' (16, 16; 0.05) (blue symbols) fore G 
{0.03, 0.04, 0.05, 0.055, 0.06, 0.08} at different orders n + 1 = 1, . . . , 12. The left-most red symbols 
are the extrapolated values. 
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FIG. 5. Self-interactions with replicas producing 1/L = l/(aNs) Coulomb terms. 
For large N$, we write 



/if \ fn(N S ) 1 



(68) 



The 0(1/ Ns) correction originates from interactions with mirror images at distances aNs, 
y/2aN s , \Z3aN s , 2aN s , . . ., see also Ref. js4]. This effectively produces a static potential 
between charges separated at distances aNs, but without self-energies (the self-energies are 
included in 5m). As illustrated in Fig. [51 the scale of such interactions is of order aNs and 
one may writd_^ 



5m(N s ) =* m --Lf3 Ua^aNg)- 1 ) + O (J^J . 
Therefore, the coefficient f n (Ns) is a polynomial of ln(A^): 

n 

f n (N s ) = J2rt ) ^(N s ), 



(69) 



(70) 



t=0 



where fn^ = f n and the coefficients fn for i > are entirely determined by f m with m < n 
and f3j (see Eq. f|44|) ). with j < n — 1. For instance, 



A) 



fi(Ns) = fi + fo^HNs, 

^Jio r /03— 2 

Z7T 87T 



f2(N S ) = f2 + 



HNs) + fo ( § ) ln 2 (iV 5j 



(71) 
(72) 



and so on. 
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There are some qualifications to this statement that we will detail below. 
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Starting at 0(a 4 ) 1 one may expect additional 0(1/ Ns) finite size terms. These would 
arise from infrared singularities of certain types of diagrams. The source of these infrared sin- 
gularities is similar to the one that results in infrared divergences of the static potential [53[ 
starting at 0(a A ). In that case these are due to the static triplet and antitriplet sources, 
which can arrange themselves into a singlet or into an octet representation at short distances 
in the pNRQCD |54( multipole expansion, giving rise to terms like a A \a.(aN s Hm) / (aN s ). In 
our case, pairs of triplet static sources can be arranged into antitriplet and sextet representa- 
tions with their mixing mediated through gluons. At higher orders different representations 
can arise, as several mirror images will interact. Nevertheless, in the limit Nt — > oo in a 
finite spatial volume, we do not expect these ultrasoft logarithms to show up. The reason is 
that aNs, besides being the typical momentum transfer between the mirror images, is also 
the infrared cut-off of the gluon momenta so that only logarithmic terms ~ \n(Ns / Ns) can 
appear. Indeed we do not detect any indication of these logs in our numerical data. 

Eq. ( |69|) can be interpreted in terms of renormalons. The fact that we can link the 
(9(l/iVs)-terin to a static potential leads to the expectation that the high order behaviors 
of f n and dominated by one and the very same renormalon. This can, e.g., be 

illustrated considering the leading dressed gluon propagator D(k) oc 1/k 2 , where k± = 0. 
With the (formal) ultraviolet cut-off 1/a and an infrared cut-off l/(aNs) this can be written 
as (ignoring lattice corrections), 



Ptxf' dkk 2 D(k) ~ - J"c n a n+1 (a- 1 ) - -L V c„a n+1 ((aA^ 1 ) , (73) 

Jl/(aN s ) a „ aN ^n 



after perturbatively expanding D(k). When re-expressing a((aNs) _1 ) in terms of a(a~ l ) we 
may consider two situations: 

(a) Ns > e n . In this limit the last term of Eq. f J73|) is exponentially suppressed in n and the 
renormalon can directly be obtained from a large order expansion of aP in powers of a. 

(b) Ns < e n . The last term of Eq. f J73|) is important and the renormalon cancels order- by- 
order in n. It is easy to visualize the importance of this term in the large-/3o limit (see the 
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discussion in Ref. 55] ). In this limit one obtains 



a 



1 \ -j2 lnn ( Ns "> ( da 



aNq J ^— ' n! \<ilna 

n>0 



large /3o 

n>0 

and therefore the large-n behaviors of c n and fn (N$) read 



(74) 



This results in the logarithms of Eq. f JTOj) to exponentiate and to cancel the 1/N$ suppression. 
Therefore, at large n, the f n /Ns terms become numerically as important as the c n -terms so 
that 

pin TVs 



5 

oo 



«=?i+i / 



■^(S"^^*' (7e) 

° v 7 i=n+l 

Note that the renormalon (n! behavior) actually cancels in the difference and that the infinite 
sum above is a convergent series (N$ minus a finite sum). 

In present-day numerical simulations, including ours, N$ < e n , and the term f n {Ns) /N$ 
needs to be taken into account, in combination with c n . A similar phenomenon was numer- 



ically observed for the static singlet energy E(r) = 2m + V(r) [50l. l56|. This teaches us that 
to correctly identify the renormalon structure of 5m, it is compulsory to incorporate the 
l/N s corrections. So far, in studies of high order perturbative expansions of the plaquette 
the corresponding finite size terms have been neglected. As we will see, our fits indeed yield 
f n ~ c n for large n, in clear support of the renormalon dominance picture. 

In the Nt — > oo limit, and up to 0{1/Ng) effects, the fit function for c n (Ns) depends on 
c n , /„ and the (3- function coefficients with i < n. In the lattice scheme only /3 , 0i and ^2 
are known. The effects of higher /3j start at 0(a 5 ). One may try to fit these, together with 
the c n and f n but their contribution cannot be resolved by the present precision of the data. 
This produces some uncertainty in our parametrization that we will add to the error. The 
reason we can neglect higher f3j in a controlled way is because the associated uncertainty 
quickly becomes negligible at high orders, once the behavior of the coefficients f n starts to 
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be governed by the d = 1 renormalon. This statement can be quantified, since we know the 
large- n behavior of the fn. Let us first consider the large-^o limit. Assuming renormalon 
dominance for the coefficients f n , we would have (see Eqs. ( 1751) and (1761) ) 



/ i— n 



UN s )=N m (^) n\J27J l ANs). (77) 

Note that terms containing higher powers of ln^s) (with % < n) are suppressed by factors 
~ 1/n! and can be neglected at large n. Therefore, for large n, the n\ factorial overcomes 
the large logs. This suppression also holds beyond the large- /3q limit o However, we may 
worry about terms with small powers of hi(Ns) (i ~ 1 in Eq. (1771)). In this case there is no 
factorial suppression, but the inclusion of fi% and /?2 can still be done in a controlled way. 
Including the running associated to 0i produces 1/n suppressed corrections to Eq. ( 1751) while 
$2 results in 1/n 2 suppression and so on. For instance in the case of f n we have 



f m = N /M T(n + i + ») 

m \2w) r(i+6) 



1 + 7 777 : 7T + O (I n 2 ) 

(n + b) (n - 1 + b) w ' 



(n 



To check this assumption and to justify the truncation at 2 we have performed separate fits 
including /3j for j < 0,1 and 2 (see Sec. IVIII) . 

Finally note that the validity of the above discussion is unaffected by any renormalon (or 
other singularity) related to the /3-function coefficients, as this would correspond to a higher 
dimension, i.e., u > 1/2 (if it existed at all). 



VI. 1/N T - AND SUBLEADING l/7V 5 -CORRECTIONS 

While most of our geometries satisfy iV^ > iV^, the dependence may still be sizable 
and cannot completely be neglected a priori. This may necessitate a combined expansion in 
powers of l/iVg and 1/Nt- The leading order (LO) correction in 1/Ng has been discussed 
in the previous section. Incorporating finite 1/Nt effects does not affect the renormalon 
structure nor the main conclusions of that section. The only subtlety that we need to revisit 
are the ultrasoft effects. In principle, a dependence on N T may appear, starting at 0(a A ). 
Nevertheless, only in the limit N T N s do we expect large logs of the type \ia(N s /N T ), as 

11 In any case, for any given n > 3, the coefficients / 7 1™' ) ^" 2 ' are completely determined by the renor- 

malization group analysis and the coefficients fx, f% and f^. 
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FIG. 6. cJ, 3 ' 0) (4,7Vt) as a function of 1/Nt compared to a linear fit. The linear term is clearly 
zero within errors. The fit gives Cq 3,0 ^(4) = 1.9221(20), to be compared to 1.92253 from DLPT, 
Eqs. 123-©. 



in this limit l/(aiVr) may act as the infrared regulator. Still our geometries are far off this 
limit. One may consider an interpolating phenomenological function like \n[N $ / {N? + Ns)\ 
between the Nt <C N$ and N$ <C Nt limits, yet the data do not seem to require these terms. 
We also stress that these terms are subleading from the renormalon point of view (d — 3). 
Therefore, in our final fit function we will not introduce them. 

We now study possible power-suppressed 1/N T effects. First we consider the low orders in 
perturbation theory. At 0(a) the fit function with finite (but large) Nt can be obtained with 
DLPT. No dependence on Nt is found. This is also confirmed by our explicit computation 
of Co with NSPT. We illustrate this for Ng = 4 in Fig. O A similar picture applies to the 
other values of N$. The leading terms in 1/N$ can also be determined using DLPT. Writing 



Am 

c {Ms, ^t) - c {l\ s > - Co 

JVC 



v m ( 1 

-^ + ( — 



(79) 
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we obtain for unsmeared coefficients and TBCxyz boundary conditions 

c o 3 dlpt = 2.1172743570834807985970 . . . , (80) 
f$°l PT = 0.76962563284(2) , (81) 
«g&r = 0-14932(3) . (82) 

Note that DLPT predicts the absence of an 0(1/ Ng) term at this order. The above result 
also applies to the adjoint source, substituting c^' p \N s ,N T ) = C A /C F cf p) (N s , N T ), where 
Ca/Cf = 2iV 2 /(iV c 2 — 1) = 9/4. We remark that f Q and v depend on the boundary 
conditions, whereas Co does not. 

As we already mentioned, the finite volume Cq'°'(Ns) = % (Ns, Nt) depend on the 



boundary conditions. It has previously been computed with PBC, original 



V in Ref. [57 1 



where intermediate semi-analytic expressions can be found, and in Refs. [33|, [34J where also 
TBCxyz and TBCxy boundary conditions were analyzed. No time dependence was found in 
either case. This absence of a time dependence at 0(a) fits with the spectral picture. The 
infinite volume coefficient was most precisely computed in Ref. (58]. Our determination of 
Co agrees with the previous results. 

At 0(a 2 ) we start to encounter a dependence on iV^. DLPT also gives us information on 
the coefficient Ci(N$, N?)- In this case we have only computed the infinite volume limit for 



59| in DLPT: 



the unsmeared coefficient using the code of Ref. 

cKlpt = C F /C A cf£l PT = 11.1425(25) . (83) 



c[ 3,0 ' ) has been computed previously in a less controlled way using finite Wilson loops 
resulting in the value cf' ^ = 11.152. In Ref. j^j agreement with this value was reported. 
Beyond 0(a 2 ) there exist no DLPT results. 

Next we address NSPT data for n > 1. We wish to understand the Nt dependence 
for Nt > Ng. For this analysis the simulations to 0(a A ) at Ns = 4 up to the very high 
Nt = 24 turn out to be particularly useful. The results are shown in Fig. [7J Note that 
in these cases the error bars are dominated by the finite Langevin timestep systematics. 
The n = 1, 2, 3 results all show the same qualitative behavior. For Nt > 10 the data are 
constant within errors, and linear fits result in slopes that are compatible with zero within 
two standard deviations. For Nt smaller than 10 we start to see a bending in 1/Nt, which 
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FIG. 7. c[ 3 2°3(4, Nt) as a function ofl/Nx, in comparison to a constant plus linear fit, a constant 
plus cubic fit, and a constant fitted only to the Nt > 10 points. 



we parameterize by a 1/Nj, function. Large powers of d are favored by the fit. The specific 
power is difficult to determine. We find a 1/iVf. fit to best describe the data, though only 
marginally better than a l/iV|> fit. Linear 1/Nt fits, however, are unsatisfactory as we can 
see in Fig. [71 We can compare the c^2°](4, oo)-values obtained averaging > 10 data vs. 
performing l/iV| fits: 9.142 ± 0.012 vs. 9.147 ± 0.010, 63.08 ± 0.13 vs. 63.16 ± 0.10 and 
508.3 ± 1.5 vs. 509.7 ± 1.0. Indeed, within our present accuracy, the large-iVr data are in 
agreement with the extrapolation. The same also holds for the smeared and octet data sets. 

We found 1/N T and 1/N T fits to also work well for other N$- values (though in these 
cases we have less data points and therefore less conclusive results). Irrespectively of the 
power d, we observe the coefficient of the 1/iV^-term to increase roughly linearly with Ng. 
From this phenomeno logical analysis we conclude that the 1/Njn effects effectively count 
as Ns/N T ~ 1/N T . Moreover, we find the coefficients of these terms to be numerically 
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FIG. 8. c± (Ng, Nt) from DLPT obtained on volumes with PBC. The fitted curves are constant 
plus cubic (1/N T ). 



small. In order to confirm this phenomenological counting, we have also explored the region 
Nt < N,s/2 for n — 1 using the PBC DLPT formulae of Ref. [stJ that only apply to 
Nt < Ng. In this case for very large volumes (N$ > 32), we indeed found a (a + bNg)/N T 
parametrization to work well, see Fig. El This means that the 1/Nt effects are clearly 
subleading, compared to the 1/Ns effects that we incorporate in our fit and also subleading 
relative to the unknown 1/Ng effects starting at O (a 5 ), due to ^3. 

While these analyses strongly indicate that the 1/Nt effects decay rapidly with Nt, the 
specific functional form is not exactly known. Therefore, our analysis strategy will be to take 
Nt > max(A^5, 11) so that the 1/Nt effects can safely be neglected. In this way we loose data 
and statistics but avoid any bias from assuming a particular functional form. To estimate 
the cut-off systematics we then vary Nt and also consider different trial fit functions. We 
discuss this issue further in the next section. 

Subleading effects of 0(1/ Ng) would be obscured by the unknown (logarithmically modu- 
lated) 1/Ns effects from higher /3-function coefficients. Therefore we will not consider these. 
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We conclude with a discussion of lattice artifacts. Formally, we may introduce an 
anisotropy a t ^ a s . In this case the lattice action, that is invariant under time or par- 
ity reversal, agrees with the continuum action up to 0(a^, a^-terms. The temporal and 
spatial lattice extents in physical units are given by a t Nx and a s Ns, respectively, so that 
the only dimensionless combinations consistent with the leading order lattice artifacts are 
af/(atNT) 2 = 1/iVy and l/Aj. Therefore, within perturbation theory, where we cannot 
dynamically generate additional scales, the LO lattice artifacts are indistinguishable from 
0(1/Nt, 1/-WJ) finite size effects which, as discussed above, are beyond our present level of 
precision. 



VII. SIMULATION RESULTS 



In this section we obtain the infinite volume coefficients of the expansions of four different 
self-energies: for fundamental and adjoint sources and using static actions with smeared 
and unsmeared time derivatives. We compare their large order behavior with theoretical 
expectations, and determine the leading renormalon normalizations N m and N m ,. We then 
convert the results into the MS scheme, using different methods, and estimate P 1 ^. 



A. Infinite volume coefficients 

Below we determine the infinite volume coefficients c^'^ defined in Eq. (|42p . Our default 
fit function for Cn (iVg, N?) (see Eq. (IMjl ) is defined in Eqs. (l68l) -(J70l). and depends on the 
fit parameters dh 1 ^ and fj R ' p ^ with j < n. This last dependence introduces a correlation 
between different n-valued coefficients, which we take into account by simultaneously fitting^ 
all c n (N s , N T ) to data up to a given order C(a nmax+1 ). By default n max + l = 20. As a sanity 
check, we have also fitted each order n independently with two fit parameters c n and /„, 
keeping the /^-values that were obtained at previous orders j < n fixed. Since this iterative 
method does not take account of all correlations, the resulting statistical errors and x 2 - y alues 
are not reliable. Nevertheless, these fits yield similar central values, illustrating that the low 

12 All the global fits to the c n (Ns,NT) data have been double checked by two different program implemen- 
tations using both Maple and Mathematica. 



34 



order coefficients are only mildly affected by higher order data. In the following we will only 
use the results of the global fits. 

To ensure that 1/Nt effects can be neglected we restrict our fits to Nt > max(Ns, vt) 
with vt = 11. In addition we restrict N$ > v$, varying v$ to explore the validity range 
of Eq. (I68p . Our "thermometer" for this will be to obtain acceptable xV-^Vdf- values and 
good agreement with Cj 3 and c^ 8 '°' ) from DLPT, Eq. ( 1831) . We find that including small 
volumes improves the quality of the fits: the values of c[ R ' ^ tend towards the expectations, 
and x 2 /Ndf, as well as the errors, are reduced. We illustrate this behavior in Table IVl We 
have observed the same behavior for different values of vt around 11, and also for the octet 
and/or smeared perturbative series. Therefore, our default setting will be vs = 4. 

TABLE V. x 2 /Nt>f, cf' 0) and cf 9 * 0) for different values of v s < N s (N T > max(N s ,ll)). The 
n = values were fixed to the DLPT result. Otherwise the x 2 / N^-p -values come out even smaller: 
1.570, 1.322, 1.209 and 1.152 respectively, whereas the coefficients barely change. The DLPT 



expectation 


is cp'°^ 


= 11.1425(25). 








vs 




9 


7 


6 


4 


X 2 /N DF 

J3,0) 




1.701 
11.120(33) 
3.919(73) 


1.431 
11.124(25) 
3.995(55) 


1.309 
11.122(17) 
4.108(36) 


1.263 
11.136(11) 
4.118(36) 



The leading parametrical uncertainty stems from the unknown 1/N$ effects associated 
to higher order terms in the /3-function: 3 , /3 4 etc., which will start affecting the fit at 
orders n + 1 > 5. As long as all singularities of the lattice /3-function in the Borel plane 
are further away than u = d/2 = 1/2 from the origin (which is the case), these higher 
Pi coefficients will not affect the leading renormalon behavior. Nevertheless, there can be 
an impact at intermediate orders. To study this, for each Cn (Ns, Nt) we perform three 
different fits, setting fix — fl 2 — (/?o), setting only f3 2 = (po,i), and using all the known 
coefficients (A), 1,2)- The resulting Cn' ^ are displayed in Table ED of the Appendix. The 
results between the (3q and /5 ,i fits start to deviate from each other significantly at n = 4 
while 0o,i,2 becomes statistically distinguishable from /3 ,i starting around n — 9. At n = 19 
there is a 25 % variance between the (3q and /3o,i,2-fits. The convergence pattern is sign 
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FIG. 9. c { „ ' >(N s )/c { r, ' ' — 1 for n G {0,1,2,3,4,5,7,9,11,15} (top to bottom). For each value of 
Ns we have plotted the data point with the maximum value of Nt- The error bars are invisible 
on the scale of the figure. The curves represent the global fit. For n = the DLPT prediction 
— (1/^s)/o 3 d°lpt/ c o 3 dlpt °f Eqs- (|80|) and ([8T]) is shown (straight line). 

alternating. The picture is similar for the smeared and the octet results. We will take the 
difference between the /?o,i and Po,i,2 results as an estimate of the error from subleading terms 
in the /3-expansion. This is our dominant source of systematic error, by far exceeding, e.g., 
our statistical errors. We remark that switching off the running altogether (/% = 0) yields a 
bad x 2 /Ndy — 3.167 (with n = fixed from DLPT) and a value of cf' ^ that is by about 20 
standard deviations away from the DLPT result. Once the running is introduced into the 
parametrization of the finite size effects, these quickly and unavoidably (see Sec. EJ grow in 
size, resulting in large cancellations with the coefficients c n . We illustrate the importance of 
this effect in Fig. [91 where we compare the fitted parametrization to the unsmeared triplet 
data on c n (Ns)/c n — 1 for various n (this also illustrates the quality of the fit). Note that 
the curvatures, i.e. the deviations from straight lines, are due to the renormalization group 
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FIG. 10. Zoom of Fig.lMf' 



or n 



running of the 1/N$ coefficients. The data clearly show the expected curvature. To illustrate 
this better, we enlarge the n = 9 curve in Fig. [TUl 

Next, we estimate the error associated to the iVy-range dependence that we have not 
accounted for in our fits. Our data run over a large variety of lattice volumes with different 
Nt- values. Our cut-off Nt > 11 eliminates a significant fraction of lattice geometries. 
However, we can still benefit from these discarded volumes, as they allow us to estimate 
the systematics associated with our choice of cut-off. We follow two strategies: i) we vary 
the cut-off v T . We display v T = 9 results in Table IVHI Reducing the cut-off increases the 
X 2 /Ndf- values, since the 1/JVy curvature is not built into our parametrization. Other than 



this there is good agreement with our v? = 11 ,#0,1,2 fits o: 
TVV-dependent term into the fit function in the following wajt£ 

Cn(Ns,N T ] 



able PVTl ii) We introduce a 



fn(N S ) ^ v n (N s ) 

Cn ~ V 



S 



13 In Ref. [10| we employed a different parameterization of the A^-dependence. The fit yielded similar results 
to those found here but using two extra parameters per order. 
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and fit to all our volumes [vt = 5). We have explored different values of d and different 
parametrizations of v n (Ng). In Sec. I VI I the low n v n (Ns) coefficients were found to increase 
with Ns- Global fits also favor this behavior. Therefore, we consider two fit functions: ii.a) 
v n (Ns) I Nj, where we construct v n (Ns) in analogy to the / n (ATg)-term, using the renormal- 
ization group running of previous orders with just one new fit parameter v n = Vn at each 
order, ii.b) v n Ns/N^, assuming a linear dependence of this term on N$. We now vary d. We 
take d — 2 for the ii.a) fit, as we obtain a good xV^df- value and agreement with c^dlpt 
within one standard deviation. Varying d increases xV^df and deteriorates this agreement. 
We take d = 3 for the ii.b) fit, as it yields a good xV^df- value and also perfect agreement 
with ci 5 dlpt- d — 2 results in a difference between the fitted value of C\ and c^dlpt of several 
standard deviations, while d = 4 and d = 5 reduce the quality of the global fit in terms of 
the x 2 -values. 

The 1/Nt effects are much less constrained by theoretical arguments than the 1/Ns 
effects. This could have resulted in a substantial increase of the number of fit parameters 
necessary to obtain acceptable xV^df- values. Fortunately, the -/V T -dependence of the data is 
much smaller than the ^-dependence. We find it remarkable that, with just one additional 
parameter per order, we can accommodate the complete Nj< dependence down to Nj< = 5. 
Note that fitting without such an iVV-term to all volumes [vt — 5, v$ — 4) we obtain an 
unacceptable xV-^df = 3.923, whereas both choices (ii.a and ii.b) yield good reduced x 2 - 
values, see Table I VI II Ansatz ii.b) gives results in perfect agreement with our — 11 
strategy, while ansatz ii.a) agrees within 1.5 standard deviations. In both cases we fixed the 
n = terms from DLPT. We notice that the coefficients v n and v n are small in size and tend 
to vanish for large n, relative to the divergent c n and /„. 

In spite of this success, we opt for the more conservative strategy of discarding data with 
Nt < Ns or Nt < 11, since our 1/iVV-fit ansatze are phenomenological and not fully under- 
stood theoretically. For the errors associated to the iVV-cut, we take the differences between 
the first columns of Tables I VI I and IVII[ as this choice is completely unbiased regarding the 



functional form of the iVr-dependence. We stress that the by far most dominant systematics 
are the unknown iV^ 1 ki l (Ns) terms. Therefore, alternative estimates of the 1/Nt effect 
would only marginally affect the final errors. 

We have completed the exploration of potential sources of systematic uncertainties. The 
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other perturbative series (smeared, octet and octet smeared) were analyzed analogously, 
with similar conclusions and precision. In particular similar x 2 -values were obtained. The 
only exception was the octet case, for which we obtained a somewhat reduced precision and 
the x 2_va hies were smaller by factors of approximately two. This could be traced to some 
geometries where the individual errors turned out much bigger. This effect then propagated 
into the final data set. 

We list the final numbers for all the infinite volume coefficients c^' 1 ^ in Table IVIIII The 
central values are taken from the first column of Table IVII The quoted errors result from 
summing statistical and theoretical uncertainties in quadrature. Schematically, we have at 
each order n 



^final = <4 a t. + °} + °T , ( 85 ) 

where ap is the difference between the first and second columns of Table IVII and <jt is the 
difference between the first columns of Tables IVT1 and I VIA We find op 3> <Jt, c s tat., so that the 
dominant error comes from logarithmic N^ 1 ln^A^-corrections, due to our lack of knowledge 
of (3^ etc.. In comparison to these unknown l/iVs-terms and the l/Ntp corrections addressed 
above, 1/iVf effects are negligible. 

In Table IVTlTl we have chosen to multiply the octet coefficients by factors Cf/Ca- In this 
normalization these will agree with the triplet coefficients for n = and n = 1 but at higher 
orders in general they will differ by 1/iV^-terms. Within our uncertainties, however, we are 
unable to resolve these differences. 

Our NSPT value cf' 0) = 11.136(11) is in good agreement with the DLPT expectation 
Eq. f )83|) . c£' ^ was calculated previously by two groups. One group determined the static 
energy, singling out the residual mass of the potential using large Wilson loops. Employing 



NSPT they obtained 4 3,0) = 86.2(0.6) (1.0) |6jJ. The second group fitted a polynomial in 



a to results of non-perturbative simulations of the Polyakov loop at various large values 



of the inverse lattice coupling (3. They obtained c^' ^ = 86.6(5) |3_4j. Our result c 



.(3,0) 



86.10(13) confirms these studies, while results for n > 2 were not known previously, e.g., 
c (3,o) = 794 5 (! 6 ) 

The same analysis also yields the 1/iVg correction coefficients fn R ' p \ where we determine 
the systematic error in the same way as for the c^' p \ We display the results in Table HXl 
For large orders the perturbative expansion should be dominated by infrared physics, 
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whereas different smearings correspond to different regularizations of the high energy behav- 
ior of the Polyakov loop. Therefore, we expect the smeared and unsmeared coefficients to 
converge to the same values for large n. This is indeed the case for the coefficients c n and f n 
of both the triplet and octet representations. Actually, the differences between smeared and 
unsmeared coefficients vanish quite rapidly, around n = 6 for the c n and already at n — 1 
for the f n . Indeed, all smeared and unsmeared values of /„ are equal within errors for both 
representations. This is to be expected, as the coefficients are related to finite size effects 
and know nothing about the specific regularization prescription for the ultraviolet behavior 
of the Polyakov loop. It is tempting to consider global fits, constraining the smeared and 
unsmeared f n values to be equal, to increase the accuracy of the results. However, to avoid 
any bias we will not explore this possibility in this article. 

We now move on to determinate the infinite volume c n /c n _i-ratios. These are obtained 
from the same fits, since we have also computed the correlation matrix. Actually, we find 
strong correlations both of the statistical and systematic errors between consecutive ex- 
pansion coefficients. Due to these correlations, the infinite volume c n /c n _i-ratios can be 
determined more precisely than the coefficients themselves. The results are displayed in 
Table [X] Up to n — 11 the errors increase. For higher orders this tendency is reversed, since 
the relative impact of the /^-value (and hence also of the unknown /3-function coefficients) 
diminishes and so do the effects of finite A^-corrections. 

As a cross-check we have also determined the coefficients c n by a direct fit to the ratio 
data 

c n -fn(N s )/N s [+v n (N s )/N«] 



-(N S ,N T ) 



Cn-1 - fn-x(N s )/N s [+v n ^(N s )/N^\ 



(86) 



Cn-1 

using the /qdlpt" vame f° r ^ ne non-smeared case and /q dlrt; obtained in the previous fit, for 
the smeared case. For the central values and error estimates we proceed in the same way as 
we did before. Doing this, overall consistent results and errors for the individual coefficients 
are found (with slightly bigger x 2 - values). The only exception is the unsmeared octet case, 
where the problems of stability that we already encountered for the c n data become magnified 
in the ratios, further reducing the precision. Subsequent c„- values are statistically correlated 
and direct fits to the ratio data take these correlations into account. We obtain similar errors 
and central values as in the previous analysis. This indicates that the statistical correlations 
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of the lattice data do not significantly affect the errors of the infinite volume coefficients, 
which are dominated by the systematics. As another related cross-check, we have computed 
the infinite volume c n /c n _i ratios using the ratio data with fit parameters c n /c n _i, f n (and 
v n or v n , see the discussion after Eq. (1841). proceeding analogously as above. From this we 
obtain very similar results to those quoted in Table Kl 

Finally, we remark that at the very high orders dominated by the renormalon behavior, 
TBC cannot and do not reduce finite volume effects, relative to PBC However, the v n - and, 
at low orders, the / n -values are significantly reduced, considerably increasing the robustness 
of the c n - and /^-determinations at intermediate and large orders. The effect is twofold. 
First, the impact of different parametrizations and of the Iow-Nt cut-off value on the 1/Nt- 
extrapolation is reduced. Second, the low-order f n times j3% and higher unknown /3-function 
coefficients that contribute to the Ng 1 ln^iVs^-terms are smaller. Therefore, the uncertainty 
due to the lack of knowledge of i > 3, becomes reduced at intermediate orders (for large 
n these effects will be small anyhow due to the renormalon dominance, see the discussion 
around Eq. (1771)) . 

B. Renormalon dominance and the determination of N m and N m . 

In the following we investigate whether the large-n behavior of the four different sets of 
c n and /„ complies with the renormalon expectation and determine the triplet and octet 
normalizations N m and N m .. 

In Fig. [HI we compare the c„/(nc„_i)-ratios summarized in Table IXl to Eq. (|59|) at differ- 
ent orders in the 1/n expansion. The LO and NLO expectations are scheme independent, 
whereas the NNLO expression depends on the scheme through (3^. For n > 8 the ratios 
clearly converge to Eq. ( 1591) . and they are within the right ball park of the NNLO prediction, 
as Fig. [H] illustrates. This is so irrespectively of the representation and smearing, confirm- 
ing the existence of the renormalon at d = 1. For completeness, we also plot the NNNLO 
0(l/n 3 ) expectation, using the 0^ estimate of Eq. (11031) . 

The renormalon picture predicts that c n ~ /„ for large n. This equality is achieved with a 
high degree of accuracy from n = 9 onwards in all four cases (compare Tables [Villi and [LX|) . 
For the values of N$ we explore, the renormalon picture also predicts a strong cancellation 
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FIG. 11. The ratios c n j '(nc n _i) /or i/ie smeared and unsmeared, triplet and octet fundamental static 
self- energies, compared to the prediction Eq. (|59p for the LO, next-to-leading order (NLO), NNLO 
and NNNLO of the 1/n expansion. For clarity, the data sets are slightly shifted horizontally by 
different off-sets. 



between c n and f n {Ns)/Ns for large n. We obtain this behavior, which we show in Fig. [H] 
with an excellent fit to the data (see, for instance, Fig. [T01 which is already at an order 
where renormalon dominance has set in). 

For each representation R we have four different sequences: Cn fl '°\ c^' 1 ^ 6 ^, fi R ' ^ and 
#' 1/6) that we may use to determine the normalizations N m (R = 3) and N m - (R = 8). To 
obtain the normalizations we divide the large n expectations Eqs. (1551) and (1571) for the triplet 
and octet representations by the coefficients obtained in Tables IVIIII and HX] respectively. We 
truncate the equations at 0(l/(n + b)) precision (NNLO), since resolving the 0(l/(n + b) 2 ) 
correction term requires the knowledge of /3g att . For large n these ratios should tend to 
constants, allowing us to extract N m and N m .. This is depicted in Figs. [12] and [13] for triplet 
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FIG. 12. N m , determined via Eq. (|55p . truncated at NNLO, from the coefficients On' , c^f' 1 , /?i 3 ' ^ 
and /n 3 ' 1 ^ 6 '. T/ie horizontal band is our final result quoted in Eq. (|89p . 



and octet sources, respectively. We use the n = 19 coefficients c^' ^ and c^' l ^\ and their 



associated errors, to obtain the normalizations (recalling that N n 



N l ^(p = 0) = 19.1(15) , C F /C A N%f(p = 0) = 18.5(16) 

rlatt / 



iVr (P = 1/6) = 18.9(15) , C F /C A N^Hp = 1/6) = 18.9(15) 



*7) 



The errors are much bigger than the differences between the four possible determinations: 
with or without smearing, using C19 or using f 19 . This is not too surprising since these 
parameters are obtained from one and the same global fit to the same data and hence highly 
correlated. Moreover, the errors are dominated by the systematics of varying the subleading 
terms of the finite volume fit function. We obtain our final result by averaging the above 
central values, with errors that accommodate both the original error bars: 



N. 



latt 



19.0 ±1.6, C F /C A N% U = 18.7 ±1.8. 
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FIG. 13. iV m . = —N_\, determined via Eq. (|57p . truncated at NNLO, from the coefficients 
Crf'°\ c^f' 1 , fn and fn' 1 ^^ . The horizontal band is our final result quoted in Eq. (|89p . To 
enable comparison with Fig. \1SX we multiply N m . by Cf/Ca- 



These numbers are included as error bands into Figs. [12] and [131 respectively. The bands 
contain all values of the n > 8 coefficients, lending credibility to our normalization estimates. 
Note that, on general grounds, we would expect the ratio Nrh/N m to differ from the Casimir 
scaling factor Ca/Cf by an 0(l/A^)-term, which naively amounts to 10%, roughly the level 
of our accuracy. We discuss this issue further in the next subsection. 

As a cross-check we also estimate the normalization from the Borel transform of the static 
energy perturbative series 
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FIG. 14. JV TO) determined from the coefficients Cn'°\ n < N, using Eq. (i92j) . T/te error band 
corresponds to the result Eq. (j89j) . 



using the function 



N 



DSP(u) = 5>£V - (1 - 2 M ) 1+b i? (Af) [H(i(«)) 



(91) 



n=0 



A^ m - (1 + ci(l - 2u) + c 2 (l - 2m) 2 + •••) + (!- 2m) 1+6 (analytic term) 
a v ' 



as it was first done in Ref. |45j for the pole mass, using ideas developed in Refs. [62|, [63|. 
Dm \u) is singular but bounded at the first IR renormalon. Therefore, we can estimate N m 
from the first coefficients of the series in u, using 



(92) 



We plot the predictions for different orders N in Fig. [TU The error is propagated from the 
error of the coefficient. Within one standard deviation the result is consistent with Eq. fl89|) 
(error band), though less precise. 
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FIG. 15. -Eg. (|93h times yfriQ, for five different values of the lattice scheme coupling constant a, 
ranging from a(u) ~ 0.096 (hq = 5) to a(u) ~ 0.036 (no = 15). The error band corresponds to the 
estimate of Eq. ([62]) . where we have used the value N m = 19.0 ± 1.6 [Eq. (|89|) /. 



Finally, we show in Fig. [15] the divergent behavior of the perturbative expansion of the 
pole mass, Eq. (TJ3]). We use the fact that r n = iscn' ^ for large n and the coefficients listed 
in Table IVIIII We compute 

2n \ ( /3 a(z/) x '' 



n -a n+l (v) 



n+l. 



v) exp 



+ 



(93) 



Al att " \(3 a(v) J \ 47T 

where we truncate Eq. fT4"9"l) at two-loop order. In Fig. [TS] we plot Eq. fl9"3"|) times ^Juq (see 
Eq. ()62])) as a function of n for a w 0.096,0.072,0.057,0.044 and 0.036. These values are 
chosen so that the minimal term in the two-loop approximation of Eq. (1601) corresponds 
to n = 5, 7, 9, 12 and 15, respectively. In terms of the inverse lattice coupling parameter 
= 3/(2vra) this covers the range 4.97 < (3 < 13.32. Orders n Q = 6, 7 (/3 w 5.8,6.6) 
are typical for present-day non-perturbative lattice simulations, with inverse lattice spacings 
1.5 GeV < a -1 < 5.2 GeV while the n = 5 value is in the strong coupling regime. 
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As expected, the contributions to the sum decrease monotonously down to an order ~ 
a n ° +1 , before starting to diverge exponentially. The horizontal error band corresponds to 
the uncertainty, estimated in Eq. (l62l) . of the sum truncated at order n 

\r I 2 3 / 2 - b 7r 3 / 2 
y^L^no+i^) = \N m \ « 1.206 \N m \ , (94) 

Alatt Po 1 l 1 + 0) 



where we used the value iV m = 19.0±1.6 [Eq. (1891 ]. Using A latt ^8.2 MeV [64J, this horizon- 
tal line corresponds to about 190 MeV. The data are very consistent with expectations, the 
only difference being that at the largest coupling (lowest scale v) the order of the minimal 
term is somewhat lower than expected (no = 3, 4 instead of no = 5). We obtain very similar 
results from the smeared and the octet data. 

At smaller a, i.e. at higher u, the minimal term c nQ a n ° +1 {v) is numerically smaller than at 
lower scales. However, this is compensated for by the linear divergence of r n = uc n , resulting 
in a similar overall uncertainty. The only difference is that to achieve this accuracy, at higher 
scales one has to expand to higher orders. 



C. Conversion to the MS scheme and determination of /3g att 

The results of the infinite volume coefficients c^' p \ N m and N m . presented above have 
been obtained in the (Wilson) lattice scheme. Translating a coefficient Cn to a different 
scheme would require the knowledge of the conversion to order a n+1 . This is completely 
beyond reach. For the case of MS, the conversion 



OmsM = «iatt(^) (1 + ^i«iatt(/i) + d 2 al tt (fi) + d 3 al tt (n) + C« tt )) 



(95) 



is known to two loops with |32|, |65|, |66j d x = 5.88359144663707(1) and |3J, |42j, |43j d 2 = 
43.4073028(2). Fortunately, only d\ is needed to determine the ratio of A-parameters, and 
and A^, since (exactly!) 



N™ S m ~ = N%%_.A laU /A m , where A^ = e *>' A 



rlatt 



latt 



28.809338139488 A 



latt • 



(96) 



This yields the numerical values 



N™ s = 0.660(56) , Cf/CaN^ 



MS 



C F /C A Nf s = 0.649(62). 



(97) 
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Other combinations of interest are (see Eqs. (|56|) and (1581) ) 

JVjp= -1.32(11), = 0.14(18) . (98) 

These results can be compared to previous determinations from continuum computations in 



the MS scheme 



50. 



56 



67 



68j . The agreement is remarkably good, which is highly nontrivial 
given the factor ~ 29 between the values of N m and in both schemes, due to the big 
difference between the A^jg- and Ai att -parameters, i.e., the large value of d\. Moreover, in the 
MS scheme the normalization was determined from the first few terms of the perturbative 
series only, while in the lattice scheme n > 9 was required. As expected, the onset of the 
renormalon dominated behavior depends on the scheme. Nowadays, several diagrammatic 
continuum perturbation theory computations in heavy quark physics have reached a level 
of precision where they become sensitive to the leading renormalon. We remark that there 
has always been some doubt about the reliability of determinations of and N^ s from 
just very few orders of perturbation theory. We have now provided an entirely independent 
determination of these objects based on many orders of the expansion that can systematically 
be improved upon. Our quenched result presented here goes beyond the present state-of-the- 
art. An analogous un-quenched determination could give similarly precise values for N^ s 



and iV^ , with direct consequences to heavy quark physics, e.g., if using the RS scheme 45]. 

To further support our conclusions, we convert the c n (N s ,N T ) lattice coefficients, and 
their ratios, into the MS scheme. As we have already mentioned, we can only exactly per- 
form this conversion up to n = 2. For n > 2 the MS coefficients and ratios will depend on 
the approximation used. If the renormalon picture is correct, the large- n ratios should be 
dominated by the renormalon behavior and all "MS-like" conversions should yield similar 
results. However, coefficients and ratios at intermediate orders will depend on the approxi- 
mation used. We consider two different MS-like conversion schemes: 
(a) MS a 

aiatt(^) = «MS (A 1 ) 1 \ 1 I A A2\ 2 ( \ ' (") 



(b) MS fc 



aiatt(^) = «is(^) (! - d i a m(^) + ( 2d l ~ rf 2)afjs(^)) • (100) 
We suspect the scheme MS a to be superior, since the translation of 1/a rather than of a 
from one scheme to another generates a renormalization group-like resummation. 
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FIG. 16. The ratio t^f'°V( nc ^i) * n the lattice and MS-like schemes, compared to the prediction 
Eq. (f59|) . NNLO and NNNLO refer to the respective MS scheme expectations. 



Our statistical data analysis 52] allows for the direct evaluation of derived/secondary 
observables. The expansion of the logarithm of the Polyakov loop is the most obvious 
secondary observable and produces the coefficients c n (N s , N T ), but we can also intertwine 
the logarithm with other functions, such as the change from the lattice to a MS-like scheme. 



MS 



We do so using Eqs. (|99l) and (11001) . In addition, we employ DLPT to obtain q 
Cf/Cac| 8 ^ = —1.3147(25), whereas the first coefficient c is scheme independent. 

In Fig. dn] we show our determination of ci 3 ' V (nc^fi ) using the two MS-like conver- 
sions (a) and (b). We only display the statistical errors associated to the fit. We have 
not performed a complete error analysis, as the MS-like conversions introduce unknown sys- 
tematics. As anticipated, both MS-like schemes converge to the renormalon expectation. 
Actually, leaving aside systematic errors, it converges to the MS NNLO expectation rather 
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than the lattice one. For the more stable MS a scheme, renormalon dominance sets in already 
at orders n ~ 5, 6, significantly earlier than in the lattice scheme. 

Following the analysis of the Sec. IVH Bl we can determine the normalization of the renor- 
malon from the coefficient C19, obtaining the estimates 

iVjP° ~ 0.51 , N^ b ~ 1.53 , (101) 

to be compared to the correct value N^ s = 0.66(6) of Eq. fl97|) . Of course, due to the 
mismatch at intermediate orders, these numbers are not trustworthy. However, considering 
the O(10 22 ) size of C19, the numbers certainly are within the right ball park. As expected, 
the scheme MS a is superior, both in terms of an earlier onset of the asymptotic behavior and 
the extracted value of the normalization N^ Sa ~ N^ s . A similar picture is obtained for all 
the other sequences except for the unsmeared octet. In this latter case the data become too 
noisy to obtain stable results. 

Since we know ft^ 13 {4]]], we can go one order higher in 1/n in the prediction for the 
ratios and coefficients of the MS-like schemes. Incorporating the running to this higher 
order into the fit function produces very small shifts of the predicted ratios and coefficients. 
This confirms that introducing consecutive orders of the /3-function into the fit leads to a 
convergent parametrization of the c n coefficients and associated ratios. 

The previous fits indicate that the asymptotic behavior of the ratios is not very sensitive 
to their values at intermediate orders. However, the normalization N m is, as the value of a 
high order coefficient c n , obtained from a global fit will, through the running of the 1/Ns 
finite size effect, also depend on n-intermediate orders. This is also so if one tries to obtain 
the coefficients c n through the ratios. The reason is that these are determined from the 

relation c n = CoH™ =1 -^-, which is sensitive to the intermediate values of In spite of 

■> c j — 1 c j - 1 

these caveats the results are encouraging and perfectly compatible with expectations. 

We expect that the renormalon dominance of the static energy expansion sets in at much 
lower orders in the MS scheme than in the lattice scheme. This is supported by the con- 
sistency of our jV m -determination with continuum estimates that are based on only a few 
orders. Also the earlier onset of the asymptotics in the MS-like schemes is coherent with 
this assumption. We can turn this argument around to estimate d 3 [cf. Eq. ( 1951) ] and (3^, 
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assuming that 



C 3,MS 



¥ ms fftV r(4 + 6) ( b 6(6-1) \ 

* m f(T+6) l 1+ (3T6) Sl+ (3 + 6)(2 + 6) S2+ J' (1 ° 2) 

Using our central value cf{^ tt = 794.5, we obtain 

d 3 ~ 365 , I3f u ~ -1.7 x 10 6 . (103) 

Eq. (I102p introduces a systematic error that is difficult to estimate. Nevertheless, we have 
checked that the value of d% varies at the few per mille level when considering the uncertainties 
of N m , 

c 3iatt> or when truncating Eq. (11021) at a lower order in 1/n. This translates into 
variations at the level of a few per cent for /3g att . We have also checked that introducing 
this estimate of fi 1 ^ in our fit function of Sec. IVH Al yields a convergent pattern (in the 
number of /3-coefficients included) for c n and c n /(nc n _i). In this case c n /(nc n _i) converges 
to Eq. f )59|) with NNNLO precision. The fit produces a somewhat smaller value of N m that 
agrees within one standard deviation with the result stated in Eq. ( |89|) . 



VIII. CONCLUSIONS 

We have determined the infinite volume coefficients of the perturbative expansions of the 
self-energies of static sources in the fundamental and adjoint representations to 0(a 20 ) in 
gluodynamics. We have employed lattice regularization with the Wilson action and two dif- 
ferent discretizations of the covariant time derivative of the Polyakov loop. The computation 
was performed using NSPT. Overall, we have obtained the infinite volume coefficients of four 
different perturbative series, which we show in Table IVIIII At high orders all series display 
the factorial growth predicted by the conjectured renormalon picture based on the operator 
product expansion. This can also nicely be seen from the normalized ratios of subsequent co- 
efficients c n /(nc n -i), which converge to Eq. ( )59l) for large n, as can be read off from Table IXl 
The coefficients that govern spatial finite size effects, /„, also grow factorially, as predicted 
by the renormalon dominance picture, see Table IIXI 

Furthermore, we have determined the normalization constant of the first infrared renor- 
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malon of a heavy quark pole mass and of the gluelump mass: 



N, 



latt 



19.0 ±1.6, Cf/CaN 1 ^ = -18.7 ± 1.8, 



(104) 



in 



N. 



MS 



0.660 ± 0.056 , C F /C A Nf s = -0.649 ± 0.062 . 



(105) 



m 



We stress that the N m - value is more than ten standard deviations different from zero, proving, 
with this significance, the existence of the d — 1 renormalon in gluodynamics. We also find 
it remarkable that we can obtain a result in a continuum scheme directly (and exactly) from 
a computation in lattice regular izat ion, with no error in the conversion. 

The above numbers are in agreement, within errors, with determinations from continuum- 
like computations, but they have been obtained using completely independent methods. In 
particular, for the first time, it was possible to follow the factorial growth of the coefficients 
over many orders, from around a 9 up to a 20 , vastly increasing the credibility of the prediction. 
The results of this article can be used to predict higher order terms of the heavy quark pole 
mass, of the static singlet and hybrid potentials and of the heavy gluino pole mass (gluelump) 
expansions. Unfortunately, at present, for the latter we do not have sufficient precision to 
discriminate Casimir scaling violation effects, suppressed by 1/AT? in the number of colors. 

Our precision is mainly limited by our knowledge of the fit function, and in particular 
of (3^. We have been able to estimate its value, f3 l ^ tt ~ —1.7 x 10 6 , assuming that the 
renormalon dominance in the MS scheme sets in around 0(a 4 ). However, an independent 
precise determination would further decrease the errors of the infinite volume coefficients 
and of the normalizations N m and N m ~. Performing simulations on larger lattice volumes 
would also be desirable, to further improve the control of finite size effects. However, the 
statistical noise increases substantially with the length of the Polyakov loop N T and we find 
simulations to become unstable for asymmetries N s ^> N T . This behavior deserves further 
study. 

While the addition of a small number of quark flavors will neither affect any of the 
qualitative conclusions presented here nor the renormalon structure of the theory, a similar 
un-quenched analysis would be very important. This would provide a reliable, independent 
determination of N^ s , including the effect of light flavors, with major impact on renormalon 
analyses in heavy quark physics and, in particular, enabling more accurate determinations 
of the heavy quark masses, including that of the top quark. 
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Appendix A: Tables 



TABLE VI. Fit of c\i ' with different approximations to the (3 -function. c ' and /q ' were 
fixed to the DLPT result. Leaving these parameters free slightly decreases the \ 2 /^Ydf '-values to 
1.111,1.152,1.177 respectively, without significant changes in any of the fit parameters. 





Po,l,2 


00,1 


Po 


/ TV T 

X /Nbf 


1.263 


1.290 


1.218 


a/10 


1.1136(11) 


1.1136(11) 


l. 1136(H) 


ca/10 


8.610(13) 


8.610(13) 


8.597(13) 


c 3 /10 2 


7.945(14) 


7.951(14) 


7.914(14) 


C4/IO 6 


8.215(26) 


8.232(26) 


8.156(26) 


C5/IO 4 


9.322(40) 


9.361(40) 


9.203(40) 


c 6 /10 b 


1.1533(61) 


1.1619(61) 


1.1292(61) 


C7/IO 7 


1.5576(96) 


1.5760(96) 


1.5067(94) 


C 8 /10 S 


2.304(16) 


2.345(16) 


2.194(15) 


c 9 /10 y 


3.747(27) 


3.837(27) 


3.499(25) 


Ii r\10 

ci /10 iU 


6.702(49) 


6.913(50) 


6.121(46) 


m /in 12 


1 31finf98"l 


1 SfiTno") 

±- .KJ\J 1 \-L\J 1 


1 1 740C8Q1 

1.11 L ±\J 1 Of ) 


C12/10 13 


2.809(24) 


2.939(24) 


2.446(21) 


C13/10 14 


6.513(56) 


6.855(58) 


5.537(51) 


C14/10 16 


1.628(14) 


1.723(15) 


1.353(13) 


C15/10 17 


4.363(38) 


4.641(40) 


3.546(33) 


Cie/10 19 


1.247(11) 


1.332(11) 


0.9925(92) 


C17/10 20 


3.785(33) 


4.059(35) 


2.953(28) 


cis/10 22 


1.215(11) 


1.308(11) 


0.930(09) 


cig/10 23 


4.118(36) 


4.446(38) 


3.094(29) 



54 



TABLE VII. Determination of c n ' using the methods i), ii.a), and ii.b) explained around Eq. (|84p 
of Sec. I VII Al The first column is the result of a fit with two parameters per order [Eq. (j68j) / to the 
Nt > max(iV,s,9) geometries. The analogous Nt > max(Ns, 11) results are displayed in the first 
column of Table I VH The second and third columns are from fits of Eq. (|84p to all volumes, with 
one extra fit parameter per order: in the second column we set d = 2 and obtain v n (Ns) from the 
renormalization group running using /3o,/3i,/?2 an d results from previous orders t^-i e ^ c -- ^ n ^he 
last column we set d = 3 and v n (Ns) = v n Ns- 





vt = 9 


v n {N s )/N T 


v n N s /N T 


X 2 /N DF 


1.666 


0.940 


1.033 


d/10 


1.1133(10) 


1.11360(89) 


1.11442(89) 


Pa/10 


8.607(12) 


8.612(10) 


8.619(10) 


c 3 /10 2 


7.940(12) 


7.944(10) 


7.947(10) 


C4/IO 3 


8.201(24) 


8.233(22) 


8.231(22) 


cs/10 4 


9.305(37) 


9.361(34) 


9.340(35) 


c 6 /10 6 


1.1512(56) 


1.1606(52) 


1.1551(53) 


c 7 /10 7 


1.5549(88) 


1.5706(81) 


1.5589(83) 


c 8 /10 8 


2.301(14) 


2.328(13) 


2.305(13) 


Cg/lO 


6. 14:2(24:) 


^70 1 /no\ 
O. (91(23 J 


6. (4b{26) 


cio/10 10 


6.695(45) 


6.790(41) 


6.695(43) 


en/10 12 


1.3144(89) 


1.3341(82) 


1.3137(85) 


C12/IO 13 


2.812(20) 


2.850(19) 


2.805(19) 


cis/10 14 


6.526(48) 


6.607(44) 


6.490(45) 


C14/IO 16 


1.632(12) 


1.652(11) 


1.620(11) 


C15/IO 17 


4.375(33) 


4.426(30) 


4.340(31) 


cie/10 19 


1.2506(94) 


1.2650(85) 


1.2401(88) 


cir/10 20 


3.796(28) 


3.839(26) 


3.764(27) 


cis/10 22 


1.2192(92) 


1.2331(83) 


1.2087(86) 


cig/10 23 


4.130(31) 


4.177(28) 


4.094(29) 
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TABLE VIII. The infinite volume 


coefficients ch' ,p " > , 


including all systematic errors. 


The unsmeared 


co-values 


are fixed using DLPT. 










cl 3 ' 0) 


c4 3,1/6) 


C ( n' 0) C F /C A 


c { t m C F /C A 


c o 


2.117274357 


0.72181(99) 


2.117274357 


0.72181(99) 


Cl 


11.136(11) 


6.385(10) 


11.140(12) 


6.387(10) 


c 2 /10 

^ / 


8.610(13) 


8.124(12) 


8.587(14) 


8.129(12) 


C3/IO 2 


7.945(16) 


7.670(13) 


7.917(20) 


7.682(15) 


C4/IO 3 


8.215(34) 


8.017(33) 


8.197(42) 


8.017(36) 


C5/IO 4 


9.322(59) 


9.160(59) 


9.295(76) 


9.139(64) 


c 6 /10 6 


1.153(11) 


1.138(11) 


1.144(13) 


1.134(12) 


C7/IO 7 


1.558(21) 


1.541(22) 


1.533(25) 


1.535(22) 


c 8 /10 8 


2.304(43) 


2.284(45) 


2.254(51) 


2.275(45) 


cg/10 9 


3.747(95) 


3.717(97) 


3.64(11) 


3.703(98) 


cio/10 10 


6.70(22) 


6.65(22) 


6.49(25) 


6.63(22) 


en/10 12 


1.316(52) 


1.306(53) 


1.269(59) 


1.303(53) 


C12/IO 13 


2.81(13) 


2.79(13) 


2.71(14) 


2.78(13) 


cis/10 14 


6.51(35) 


6.46(35) 


6.29(37) 


6.45(35) 


C14/IO 16 


1.628(96) 


1.613(97) 


1.57(10) 


1.614(97) 


C15/10 17 


4.36(28) 


4.32(28) 


4.22(29) 


4.33(28) 


cie/10 19 


1.247(86) 


1.235(86) 


1.206(89) 


1.236(86) 


cit/10 20 


3.78(28) 


3.75(28) 


3.66(28) 


3.75(28) 


cis/10 22 


1.215(93) 


1.204(94) 
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TABLE IX. The 1/Ns correction coefficients fh , including all systematic errors. The unsmeared 



fo-values are fixed using DLPT. 
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TABLE X. The infinite volume ratios c n / \nc n l\j, including all systematic errors. Note that 
A)/(2tt) « 1.7507. 
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